/**
 * Consider the modified harmonic series by omitting every term where 
 * the denominator has three or more consecutive equal digits. This 
 * series converges. 
 *
 * Find the value the series converges to. Give your answer rounded to 
 * 10 digits behind the decimal point.
 */

#include <iostream>
#include <iomanip>
#include <numeric>
#include <vector>
#include "euler/matrix.hpp"
#include "euler/binomial.hpp"
#include "euler.h"

BEGIN_PROBLEM(368, solve_problem_368)
	PROBLEM_TITLE("A Kempner-like series")
	PROBLEM_ANSWER("253.6135092068")
	PROBLEM_DIFFICULTY(3)
	PROBLEM_FUN_LEVEL(3)
	PROBLEM_TIME_COMPLEXITY("")
	PROBLEM_SPACE_COMPLEXITY("")
END_PROBLEM()

template <class Ty, int Base, int N1>
Ty kempner(const int (&T)[N1][Base], Ty eps, bool verbose = false)
{
	const int N = N1 - 1;

	typedef euler::matrix<Ty,N,1,1,1> VecN; // one-based column vector.
	typedef euler::matrix<Ty,N,N,1,1> MatN; // one-based square matrix.

	auto sum = euler::msum<VecN>;
	auto a = [](int k, int w, int m) -> Ty {
		Ty c = std::pow(-1.0, (double)w) * euler::binom<Ty>(k + w - 1, w);
		return std::pow((Ty)m / Base, w) * c / std::pow((Ty)Base, k);
	};

	if (verbose)
	{
		std::cout << "eps = 1e" << (int)std::log10(eps) << std::endl;
	}

	// Compute the extrapolation matrix A and B.
	// A(j,l) = number of digits which, when appended to class l, lead to class j.
	// B = inv(I-A) - I.
	// b = sum(B) by column, then transposed.
	MatN A = {0};
	for (int l = 1; l <= N; l++)
	{
		for (int m = 0; m < Base; m++)
			if (T[l][m])
				++A(T[l][m], l);
	}
	A /= Base;
	MatN I = MatN::identity();
	MatN B = euler::inv(I-A) - I;
	VecN b = {0};
	for (int l = 1; l <= N; l++)
		for (int j = 1; j <= N; j++)
			b(l) += B(j, l);

	// Compute the sum of 1/s for single digit denominators.
	Ty total = 0;
	for (int m = 1; m < Base; m++)
	{
		if (T[0][m])
		{
			Ty s = (Ty)m;
			total += 1 / s;
		}
	}
	if (verbose)
	{
		std::cout << "1 digits: Sum=" << total << std::endl;
	}

	// Compute the sum of 1/s^k for double-digit denominators.
    // Compute for each k until the sum of 1/s^k < eps.
	std::vector<VecN> psi;
	for (int k = 1; ; k++)
	{
		VecN t = {0};
		for (int m1 = 1; m1 < Base; ++m1)
		{
			for (int m2 = 0; m2 < Base; ++m2)
			{
				if (T[0][m1] && T[T[0][m1]][m2])
				{
					Ty s = Ty(m1 * Base + m2);
					int j = T[T[0][m1]][m2];
					t(j) += std::pow(s, -k);
				}
			}
		}
		if (k > 1 && sum(t) < eps)
            break;
		psi.push_back(t);
	}
	total += sum(psi[0]);
	if (verbose)
	{
		std::cout << "2 digits: Sum=" << sum(psi[0]) 
			<< ", K=" << psi.size() << std::endl;
	}

	// Compute psi_i^j(k) iteratively until k = 1.
	for (int i = 3; psi.size() > 1; ++i)
	{
        // Compute each psi_i^j(k).
		for (int k = 1; k <= (int)psi.size(); k++)
		{
			VecN t = {0};
			for (int l = 1; l <= N; l++)
			{
				for (int m = 0; m < Base; m++)
				{
					if (T[l][m])
					{
						int j = T[l][m];
						for (int w = 0; w <= (int)psi.size() - k; w++)
							t(j) += a(k,w,m)*psi[k-1+w](l);
					}
				}
			}
			if (verbose)
			{
				// std::cout << "Psi_" << i << "(" << k << ") = " << sum(t) << std::endl;
			}
			if (k > 1 && sum(t) < eps)
			{
				psi.resize(k-1);
				break;
			}
			else
			{
				psi[k-1] = t; // overwrite psi in-place.
			}
		}

		// Debug output.
		if (verbose)
		{
			std::cout << i << " digits: Sum=" << sum(psi[0]) 
				<< ", K=" << psi.size() << std::endl;
		}

		// Accumulate the sum of i-digit denominators.
		total += sum(psi[0]);
	}

	// Extrapolate the sum.
    Ty remaining = euler::inner_product<Ty>(b, psi[0]);
	total += remaining;
	return total;
}

template <class Ty, int Base, int N1>
Ty kempner(const int (&T)[N1][Base], int accuracy, bool verbose = false)
{
	// Use eps = 1/10^(accuracy+2+N1).
	Ty eps = 1/std::pow((Ty)Base, accuracy+2+N1);
	return kempner(T, eps, verbose);
}

static void solve_problem_368()
{
	const int Accuracy = 10;
	std::cout.setf(std::ios::fixed, std::ios::floatfield);
	std::cout.precision(Accuracy);
	bool verbose = false;

#if 0 // Tests
	
	// Not containing 9.
	if (1)
	{
		int T[2][10] = {
			{ 0, 1, 1, 1, 1, 1, 1, 1, 1, 0 }, // single digit class
			{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 0 }, // not ending with 9
		};	
		std::cout << "Result: " << kempner<double>(T, Accuracy, verbose) << std::endl;
		std::cout << "Verify: 22.920676619264150" << std::endl;
	}

	// No odd digits.
    if (1)
	{
		int T[2][10] = {
			{ 0, 0, 1, 0, 1, 0, 1, 0, 1, 0 }, // single digit class
			{ 1, 0, 1, 0, 1, 0, 1, 0, 1, 0 }, // containing no odd digits
		};
		std::cout << "Result: " << kempner<double>(T, Accuracy, verbose) << std::endl;
		std::cout << "Verify: 1.962608412994616" << std::endl;
	}

	// Not containing 42.
    if (1)
	{
		int T[3][10] = {
			{ 0, 1, 1, 1, 2, 1, 1, 1, 1, 1 }, // single digit class
			{ 1, 1, 1, 1, 2, 1, 1, 1, 1, 1 }, // not ending with 4 (nor 42)
			{ 1, 1, 0, 1, 2, 1, 1, 1, 1, 1 }, // ending with 4
		};
		std::cout << "Result: " << kempner<double>(T, Accuracy, verbose) << std::endl;
		std::cout << "Verify: 228.446304159230813" << std::endl;
	}

#else

	int T[21][10] = {
		{ 0, 2, 3, 4, 5, 6, 7, 8, 9, 10 }, // single digit class
        {11, 2, 3, 4, 5, 6, 7, 8, 9, 10 }, // ending with 0 but not 00
        { 1,12, 3, 4, 5, 6, 7, 8, 9, 10 }, // ending with 1 but not 11
        { 1, 2,13, 4, 5, 6, 7, 8, 9, 10 }, // ending with 2 but not 22
        { 1, 2, 3,14, 5, 6, 7, 8, 9, 10 }, // ending with 3 but not 33
        { 1, 2, 3, 4,15, 6, 7, 8, 9, 10 }, // ending with 4 but not 44
        { 1, 2, 3, 4, 5,16, 7, 8, 9, 10 }, // ending with 5 but not 55
        { 1, 2, 3, 4, 5, 6,17, 8, 9, 10 }, // ending with 6 but not 66
        { 1, 2, 3, 4, 5, 6, 7,18, 9, 10 }, // ending with 7 but not 77
        { 1, 2, 3, 4, 5, 6, 7, 8,19, 10 }, // ending with 8 but not 88
        { 1, 2, 3, 4, 5, 6, 7, 8, 9, 20 }, // ending with 9 but not 99
        { 0, 2, 3, 4, 5, 6, 7, 8, 9, 10 }, // ending with 00
        { 1, 0, 3, 4, 5, 6, 7, 8, 9, 10 }, // ending with 11
        { 1, 2, 0, 4, 5, 6, 7, 8, 9, 10 }, // ending with 22
        { 1, 2, 3, 0, 5, 6, 7, 8, 9, 10 }, // ending with 33
        { 1, 2, 3, 4, 0, 6, 7, 8, 9, 10 }, // ending with 44
        { 1, 2, 3, 4, 5, 0, 7, 8, 9, 10 }, // ending with 55
        { 1, 2, 3, 4, 5, 6, 0, 8, 9, 10 }, // ending with 66
        { 1, 2, 3, 4, 5, 6, 7, 0, 9, 10 }, // ending with 77
        { 1, 2, 3, 4, 5, 6, 7, 8, 0, 10 }, // ending with 88
        { 1, 2, 3, 4, 5, 6, 7, 8, 9, 0  }, // ending with 99
	};

#if 0
	std::cout << "Result: " << kempner<double>(T, Accuracy, verbose) << std::endl;
	if (verbose)
		std::cout << "Verify: 253.6135092068" << std::endl;
#else
	std::cout << kempner<double>(T, Accuracy, verbose) << std::endl;
#endif

#endif
}
